• chemical modification to an individual DNA base (Cytosine)
• doesn’t change coding sequence

image source: https://www.epigentek.com

• heritable upon cell division
• dynamic

image source: https://slideplayer.com/slide/10719924/

Laird 2010 (https://doi.org/10.1038/nrg2732)

• Bisulfite sequencing (e.g. WGBS, RRBS)
  • basepair resolution
  • binomial count of methylated/unmethylated reads
• Analysis: beta-binomial regression
• Packages: bsseq, dss, dmrseq

• QC (fastQC/multiQC)
• adapter/quality trimming Trim Galore!
• alignment with Bismark
• read into R with bsseq package

• differentially methylated cytosine (DMC) or locus (DML):
  • array: test each probe
  • bisulfite sequencing: test each cytosine
• differentially methylated region (DMR):
  • array: test groups of neighboring probes
  • bisulfite sequencing: test groups of cytosines

• For cytosine \(i\) and sample \(j\) in condition \(s\), we have:
• \(M_{ij}\) reads corresponding to methylation
• \(C_{ij}\) total reads covering \(i\)
• \(p_s\) be the methylation probability in condition \(s\)
• Let \(M_{ij} | C_{ij} \sim Binom(C_{ij}, p_{s})\)
• How to test whether \(p_1 = p_2\)?

• Generalized linear model for probability of success \(p\): \[log(\frac{p}{1-p}) = \boldsymbol{X\beta}\]
• Link function \(g(p) = log(\frac{p}{1-p})\) describes relationship between proportional response and linear predictor
• No closed form; fit with iterative ML estimation
• Interpret coefficients on original scale with inverse link function \(g^{-1}\): \[ p = \frac{e^{\boldsymbol{X\beta}}}{1+e^{\boldsymbol{X\beta}}} \]
• In R: glm(cbind(successes, failures) ~ x, family="binomial")

• In ordinary logistic regression, \(p\) is assumed to be constant for all samples with in group
• To model overdispersion, we might want to allow this quantity to vary
• For example, let \[p_s \sim Beta(\alpha_s, \beta_s)\]
• Overdispersion parameter \(\phi = \frac{1}{\alpha_s + \beta_s + 1}\) contributes to increased variance
• In R: aod::betabin(cbind(successes,failures) ~ x, ~ x)

• What happens to link function \(log(\frac{p}{1-p})\) if \(p=0\) or 1?
  • binomial regression unstable for fully methylated or unmethylated cytosines
• Computationally intensive to fit model at every cytosine
• DSS: Park & Wu 2016 (https://doi.org/10.1093/bioinformatics/btw026)
  • Differential methylation under general experimental design
  • Alternate link function: \(arcsine(2p-1)\)
  • Approximate fitting with Generalized Least Squares (GLS)

• Hybrid of linear regression and generalized linear regression
• Pro: stable & closed form estimates (fast)
• Con: approximate
• Key idea: flexible covariance structure allows for specification of approximate beta-binomial error

Examples: BSmooth (Hansen et al. 2012), DSS (Park et al. 2016)

\[ \text{False Discovery Rate (FDR)} = E\Big[\frac{FP}{FP + TP}\Big]\]

• \(FDR_{CpG} = 2/10 = 0.20\)
• \(FDR_{DMR} = 1/2 = 0.50\)

• Korthauer et al. 2018 (https://doi.org/10.1093/biostatistics/kxy007)
• Key ideas:
  • model methylation signal over region
  • permutation to acheive accurate FDR control
• dmrseq: accurate inference for detection of differentially methylated regions (Bioconductor)

Korthauer et al. 2018 (https://doi.org/10.1093/biostatistics/kxy007)

• “Promoter DNA methylation is generally not sufficient for transcriptional inactivation”

• original study used DSS (approach grouping individual CpGs)
• main question : does accurate inference of methylation increase make a difference in the conclusion?

Korther & Irizarry 2018 (https://doi.org/10.1101/381145)

Korther & Irizarry 2018 (https://doi.org/10.1101/381145)

Korther & Irizarry 2018 (https://doi.org/10.1101/381145)